TP- Démonstration du filtrage et de la manipulation du Spectrogramme d'un signal de parole¶

I. Fonctions utiles pour le calcul du spectrogramme et le filtrage¶

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import Audio, display
import warnings
warnings.filterwarnings("ignore")


#Fenêtre de hamming
def hamming(T): return 0.54-0.46*np.cos(2*np.pi*np.arange(T)/(T-1))

#Calcul du spectre amplitude et phase du signal data, taille de la fenetre=T, pas=p
#mettre pre=True dans l'appel pour activer la préaccentuation
#mettre ham=True pour activer le fenêtrage par hamming
#Tfft=taille de la fft, si > T le zéro padding sera activé(voir cours acoustique)
def spectrogram(data,fs,T=512,p=32,pre=False,ham=True,Tfft=None,norm=True):
    if Tfft is None: Tfft=T #si la taille de la fft n'est pas spécifiée prendre Tfft=T        
    if norm: data=(data-np.mean(data))/np.std(data) #normaliser le signal sur son écart type
    
    if pre: data[1:]-0.97*data[:-1]      #préaccentuation
    s=[data[i:i+T] for i in range(0,len(data)-T,p)] # fenêtrage
    if ham : s=s*hamming(T)          # multiplication par hamming            
    s=np.fft.fft(s,Tfft)                 #Transformée de Fourier
    s=s[:,:int(Tfft/2)]                  #couper le spectre en 2 pour éliminer l'effet mirroir
    
    #retourner les spectres d'amplitude et de phase
    return {"ampl": np.abs(s), "phase": np.angle(s), "fs":fs, "duree":len(data)/fs, "T": T, "p": p, "Tfft": Tfft}


#Cette fonction affiche un spectrogramme d'amplitude
def afficher(spec):
    fig, ax1 = plt.subplots()
    ax2 = ax1.twinx()
    ax3 = ax1.twiny()
    freq=np.linspace(0, spec["fs"]/2, int(spec["Tfft"]/2)) #labels de l'axe des fréquences
    temps=np.linspace(0, spec["duree"], len(spec["ampl"])) #labels de l'axe du temps
    ax1.pcolormesh(temps, freq, np.log(spec["ampl"]+1).T, cmap='gray_r')  #afficher le spectre avec correction des axes
    ax1.set_ylabel('Fréquence (Hz)')
    ax1.set_xlabel('Temps (sec)')
    ax2.set_yticks(np.round(ax2.get_yticks()*int(spec["Tfft"]/2)))
    ax2.set_ylabel('Indices des colonnes dans la matrice spec')
    ax3.set_xticks(np.round(ax3.get_xticks()*len(spec["ampl"])))
    ax3.set_xlabel('Indices des lignes dans la matrice spec')
    plt.show()
    

#Retour à partir de l'amplitude et de la phase au signal temporel 
#Ne pas utiliser la préaccentuation si vous voulez utiliser cette fonction
def spec2wav(spec):
    p=spec["p"]; T=spec["T"]
    ne=(len(spec["ampl"])-1)*p+T  #estimation le nombre d'échantillons du signal
    signal=np.zeros(ne)  #initialiser le signal par des zéros 
    trams=np.zeros(ne)  #initialiser le nombre de trames par des zéros 
    
    temp=spec["ampl"]*np.exp(1j*spec["phase"]) #recombiner l'amplitude avec la phase (nombre complexe)
    temp=np.fft.ifft(temp, spec["Tfft"]) #retour au domaine temporel par une fft inverse
    temp=np.real(temp) #ne garder que la partie réelle
    
    for i in range(len(temp)): #fenêtrage inverse
        signal[i*p:i*p+T]+=temp[i,:T]
        trams[i*p:i*p+T]+=1

    return signal/trams, spec["fs"] #retourner le signal reconstitué et la fréquence d'échantillonnage

Chargement d'un fichier wav et affichage de son spectrogramme¶

In [2]:
dataset_path="Desktop\\TP dataset\\dataset\\" #le chemin vers le dataset
filename="7_george_0.wav" #le nom du fichier à charger

print('Chargement de :', filename)#chargement d'un fichier audio wav
fs, data = wavfile.read(dataset_path+filename) #fs: fréquence d'échantillonnage 

print('Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))
plt.plot(np.arange(len(data))/fs,data); plt.show() #affichage du signal temporel

spec=spectrogram(data, fs, T=256, Tfft=1025) #calcul du spectrogramme 
afficher(spec)  #affichage du spectrogramme
Chargement de : 7_george_0.wav
Freq échantillonnage: 8000 Hz, durée: 0.641375 sec
Your browser does not support the audio element.

Reconvertir le spectrogramme en signal temporel¶

In [3]:
newdata, fs= spec2wav(spec)
display(Audio(newdata,rate=fs))
Your browser does not support the audio element.

II. Filtrage¶

Filtre passe bas¶

In [4]:
spec=spectrogram(data, fs, T=256, Tfft=1024) #calcul et affichage d'un spectrogramme avec fenêtre de hamming
spec["ampl"][:,128:]=0  #annuler toutes les fréquences > 1000Hz
afficher(spec) #affichage du spectrogramme
newdata, fs= spec2wav(spec)  #retour au domaine temporel
display(Audio(newdata,rate=fs)) #écouter le résultat
Your browser does not support the audio element.

Filtre passe haut¶

In [5]:
spec=spectrogram(data, fs, T=256, Tfft=1024) #calcul et affichage d'un spectrogramme avec fenêtre de hamming
spec["ampl"][:,:256]=0  #annuler toutes les fréquences < 2000Hz
afficher(spec) #affichage du spectrogramme
newdata, fs= spec2wav(spec)  #retour au domaine temporel
display(Audio(newdata,rate=fs)) #écouter le résultat
Your browser does not support the audio element.

Filtre passe bande¶

In [6]:
spec=spectrogram(data, fs, T=256, Tfft=1024) #calcul et affichage d'un spectrogramme avec fenêtre de hamming
spec["ampl"][:,:128]=0  #annuler toutes les fréquences < 1000Hz
spec["ampl"][:,256:]=0  #annuler toutes les fréquences > 2000Hz
afficher(spec) #affichage du spectrogramme
newdata, fs= spec2wav(spec)  #retour au domaine temporel
display(Audio(newdata,rate=fs)) #écouter le résultat
Your browser does not support the audio element.

Questions¶

Q1. Essayez les filtres en haut sur les mots two, four, seven, eight et notez les phonèmes qui sont affectés par le filtrage.

Q2. Chargez le fichier mots_bruit.wav et essayez d'y filtrer le bruit avec la méthode de soustraction spectrale.

In [7]:
fs, data = wavfile.read("mots_bruit.wav") #chargement d'un fichier audio wav
print('Signal orginal: Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs)) 
Signal orginal: Freq échantillonnage: 8000 Hz, durée: 19.569875 sec
Your browser does not support the audio element.

III. Manipulation de la fréquence fondamentale (f0) par interpolation¶

In [8]:
from scipy.interpolate import interp1d

# Etirer le spectre d'amplitude avec le facteur "a"  (new f0= f0*a)  
def interpolate(spec, a):
    n,m=spec.shape
    newspec=np.zeros((n,m))
    
    for i in range(n):
        f=interp1d(np.arange(m), spec[i,:])
        new_x=np.arange(m)/a
        new_x[new_x>m-1]=m-1
        newspec[i,:]=f(new_x)
    return newspec


######## début du programme ####################
fs, data = wavfile.read("fr.wav") 
print('Signal orginal: Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))


a=0.7
spec=spectrogram(data, fs, p=16, T=512) 
spec["ampl"]=interpolate(spec["ampl"], a)
newdata, fs= spec2wav(spec)
print('Signal avec f0 manipulé de %f'%a)
display(Audio(newdata,rate=fs))
      
a=1.3
spec=spectrogram(data, fs, p=16, T=512) 
spec["ampl"]=interpolate(spec["ampl"], a)
newdata, fs= spec2wav(spec)
print('Signal avec f0 manipulé de %f'%a)
display(Audio(newdata,rate=fs))
Signal orginal: Freq échantillonnage: 16000 Hz, durée: 19.93 sec
Your browser does not support the audio element.
Signal avec f0 manipulé de 0.700000
Your browser does not support the audio element.
Signal avec f0 manipulé de 1.300000
Your browser does not support the audio element.

Questions¶

Q1. Essayez la manipulation sur les fichiers fr.wav, ar.wav, en.wav avec différents paramètres "a" et déduire l'impact de la valeur de "a" sur la parole

IV. Simulation du filtrage de l'oreille humaine (Echelle Mel)¶

IV.1 Fonctions utiles¶

In [9]:
## Fonctions utiles de conversion
def Mel2Hz(mel): return 700 * (np.power(10,mel/2595)-1)
def Hz2Mel(freq): return 2595 * np.log10(1+freq/700)
def Hz2Ind(freq,fs,Tfft): return (freq*Tfft/fs).astype(int)

#Réalisation de nf filtres sur l'échelle Mel d'une fréquence min à une fréquence max
def FiltresMel(fs, nf, Tfft, fmin, fmax):
    Indices=Hz2Ind(Mel2Hz(np.linspace(Hz2Mel(fmin), Hz2Mel(min(fmax,fs/2)), nf+2)),fs,Tfft)
    filtres=np.zeros((int(Tfft/2), nf))
    for i in range(nf):
        f=Indices[i+2]
        if (Indices[i]-f)%2==0: f=min(f+1,int(Tfft/2))
        if (f-Indices[i])>1: filtres[Indices[i]:f,i]=hamming(f-Indices[i])
        else: filtres[Indices[i]:f,i]=1
    return filtres

IV.2 Construction de filtres auditifs sur l'échelle Mel et filtrage d'un spectrogramme¶

In [10]:
nf=20    #nombre de filtres
T=256
Tfft=1024
fs=8000

#réalisation des filtres auditifs
filtres=FiltresMel(fs, nf=nf, Tfft=Tfft, fmin=500, fmax=3500)
plt.plot(filtres) 
#affichage des filtres
plt.xticks(plt.xticks()[0][1:-1],np.round(plt.xticks()[0][1:-1]*fs/Tfft))
plt.xlabel("Fréquence (Hz)")
plt.title("%d Filtres auditifs"%nf)
plt.show()

## Filtrage 
filename="7_george_0.wav"
fs, data = wavfile.read(dataset_path+filename) 
spec=spectrogram(data, fs, T=T, Tfft=Tfft) #calcul du spectrogramme
afficher(spec);  #affichage du spectrogramme

#affichage de la sortie des filtres auditifs
spec_filtre=(spec["ampl"]@filtres)/10
plt.pcolormesh(np.log(spec_filtre+1).T, cmap="gray_r")
plt.title("Spectrogramme filtré par %d filtres auditifs"%nf)
plt.ylabel("Sortie de chaque filtre")
plt.show()

IV.3 Compression de la parole par les filtres auditifs¶

Objectif: Insérer un maximum de zéros en éliminant les sons qui ne sont pas perçus par l'oreille humaine

In [ ]:
nf=36    #nombre de filtres
T=512
fs, data = wavfile.read("fr.wav") #chargement d'un fichier audio wav

#Foncrion Compression du signal par les filtres auditifs
def Compression(spec, nf, fmin=200, fmax=4000):
    filtres=FiltresMel(spec["fs"], nf, spec["Tfft"], fmin, fmax)
    newspec=np.zeros(spec["ampl"].shape)
    for i in range(filtres.shape[1]):
        ind=np.argmax(spec["ampl"]*(filtres[:,i]>0),axis=1)
        for j in range(len(ind)):  #ne grader que la valeur max de chaque bande de filtre
            newspec[j, ind[j]]=spec["ampl"][j, ind[j]]
    return newspec


################# début du programme #################
print('Signal orginal: Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))
spec=spectrogram(data, fs, p=int(T/32), T=T) #calcul du spectrogramme 
afficher(spec)
plt.show()

newspec=Compression(spec, nf=nf, fmin=0, fmax=8000)  #compression
spec["ampl"]=newspec  #remplacer le spectre d'amplitude par le spectre compressé
newdata, fs= spec2wav(spec)  #retour au domaine temporel
print("Signal compressé avec %d filtres"%nf)
display(Audio(newdata,rate=fs)) #écouter le résultat
afficher(spec)
Signal orginal: Freq échantillonnage: 16000 Hz, durée: 19.93 sec
Your browser does not support the audio element.
Signal compressé avec 36 filtres
Your browser does not support the audio element.

Questions¶

Sachant que plus le nombre de filtres est petit, plus le son est compressé.
Q1. Essayez plusieurs paramètres nf, T, fmin, fmax, sur les fichier fr.wav, ar.wav, en.wav afin d'obtenir un bon taux de compression avec une bonne qualité de son.